In ultrasonic imaging systems, Doppler effect has been widely used in blood flow detection. For example, in an ultrasonic system for examining heart, artery and vein by the use of ultrasonic Doppler technique, it is necessary to extract related parameters from the Doppler spectrum diagram to estimate the dynamic states of blood flow of the heart and vessels. Those parameters include peak velocity, diastolic minimal velocity, resistance index, systolic/diastolic index, and etc. The determination of many of those parameters depends on the peak velocity and mean velocity of blood flow. While the detection of the peak velocity and mean velocity of blood flow can be converted to the detection of the maximum frequency shift and the mean frequency shift of the ultrasound echo. The principle is that an ultrasonic transducer transmits ultrasonic waves into the body of the person to be examined, and receives the ultrasonic echoes returned from the body of the person; the frequency of the ultrasonic echo shifts due to the scattering of the ultrasonic waves by the moving red blood cells in the vessel, and the magnitude of the frequency shift is related to the central frequency of the transmitted ultrasonic wave and the moving velocity of the red blood cells, thus the moving states of the red blood cells can be obtained by calculation, so long as the frequency shifts of the ultrasonic echoes are measured.
FIG. 1 illustrates a basic flowchart of signal processing in a conventional ultrasound Doppler spectrum analysis system. The ultrasonic echo signal forms an RF echo signal after beam synthesis, and the RF echo signal is decomposed into two components by a demodulation module, one is an I (In-phase) component and the other one is a Q (Quadrate) component. In a pulsed Doppler signal processing system, the I and Q components are range-gated by the system, respectively, that is, the two component signals are accumulated respectively in specific time intervals; the operator selects the accumulation time interval and the length of the pulse transmitted by the pulsed Doppler system according to actual situations; while the range gating procedure is not required by a continuous wave Doppler signal processing system. The I and Q components (for the pulsed Doppler signal processing system, refer to the two range accumulated components) are further filtered by a wall filter (high pass filter), respectively, to remove the clutters caused by static or lowly moving tissues and to acquire the two component signals mainly including the echoes caused by the motion of red blood cells; and these two component signals are sent to a spectrum estimation module. The spectrum estimation module generally estimates power spectrum by the use of fast Fourier transform (FFT); and the number of points of the FFT may be 128 or 256 as well. Since the dynamic range of the estimated power spectrum is too broad, each of the estimated power spectrum should be compressed by the system to be fit into the gray scale display range, and the power spectral intensity of the corresponding time and velocity (frequency shift) is displayed on the Doppler spectrum diagram on the screen.
Since the detection of the peak velocity of blood flow is in correspondence with the detection of the maximum Doppler frequency shift, a conventional method of manual detection of peak velocity of blood flow by the use of ultrasonic Doppler spectrum diagram comprises the steps of: reserving the power spectrum diagram of several cardiac cycles first; manually tracing the peak values of the spectrum by the operator according to the power spectrum diagram on the display screen; and finally, calculating each of the parameters by a computer according to said traced peak values. The disadvantages of such manual detection are that the tracing of peak velocities by the operator is tedious and time consuming, the repeatability is poor, and the estimation accuracy is low; furthermore, the operator must stop the acquisition of Doppler signals in order to trace peak velocities during detection, thus it is unable to estimate in real time.
It is, therefore, necessary for the system in FIG. 1 to further include an automatic envelope detection module, such that the variations of the peak and mean velocities of blood flow versus time can be automatically traced and displayed on the Doppler spectrum diagram in real time. The automatic detection can be performed on the estimated power spectrum after the spectrum estimation, or it can be performed on the compressed power spectrum as shown in FIG. 1. Thus the peak velocity, mean velocity, and other related parameters of blood flow can be obtained automatically by the cardiologist in real time.
In principle, it seems relatively simple to detect the Doppler peak velocity of blood flow, and only detection of the maximum Doppler frequency shift is necessary. In fact, the detection is affected by two main factors: one factor is the inherent widening of the acoustic spectrum, because the number of data points used to estimate the power spectrum is limited (such as 128 or 256 points), the estimated bandwidth of frequency spectrum is wider than the ideal one, the power spectrums are distributed in the whole cut-off frequency band; the effect of this factor on the maximum Doppler frequency shift is difficult to be quantitatively determined. The another factor is the noise contained in the Doppler signal per se, a signal and a noise are often included in the Doppler spectrum diagram, a turning point between the signal and the noise can be found by detecting the maximum frequency shift; however, the turning point of the spectrum from noise to signal is not very obvious, the Doppler frequency shift varies rapidly during the blood ejection period of one cardiac cycle, but the variations during other periods are slow; in addition, for a specific frequency, the signal-to-noise ratio of the spectrum diagram varies with time. Therefore, various methods for automatically detecting peak velocity and mean velocity of a blood flow have been successively proposed in an attempt to overcome the above mentioned affections.
In order to accurately and steadily estimate the envelope of Doppler spectrum diagram, a boundary differentiating the Doppler spectral signal from the noise on the spectrum diagram is necessary. One of the differentiating methods is to set a threshold. Those on the spectrum diagram greater than the threshold are regarded as signals, and those lower than the threshold are regarded as noises. The threshold may be set to a certain fixed percentage of the sum of all signals and noises. Under the condition of high signal-to-noise ratio (SNR), the effect of this method is preferable. In clinical applications, this method is severely affected by SNR and bandwidth. In the case of low SNR, the noise threshold set at a fixed percentage may be lower than the actual noise level; while in the case of high SNR, the noise threshold set at a fixed percentage may be higher than the actual noise level. Because the noise estimation according to this method is realized by averaging the spectral lines adjacent to a cut-off frequency, therefore when the SNR is relatively low, the estimated peak value shifts in the positive direction; when the SNR is relatively high, the estimated peak value shifts in the negative direction. The mean threshold according to this method may vary with each of the spectral lines due to the randomness of noise, which results in variation of the estimated threshold level therewith. However, the actual mean noise level is relatively stable, thus the method for automatically tracking envelope may cause unpredicted results.
In an article “Comparison of four digital maximum frequency estimators for Doppler ultrasound”, Ultrasound in Med. & Biol., Vol. 14, No. 5, pp. 355-363, (1988), Larry Y. L. et al made a comparison of four methods of estimating the maximum frequency shift of Doppler frequency, and proposed an improved percentage method (referred to as “Mixed Method”). This method comprises firstly calculating a spectral integral curve for each of the spectral lines, then analyzing the integral spectral line. Generally, the spectrum diagram energy is mainly concentrated in the lower frequency portion, thus the spectral integral curve varies quickly in the potion of lower frequencies, and slowly in the portion of higher frequencies. This mixed method comprises searching for the intersection of a predetermined straight line with the spectrum integral curve, and regards the frequency corresponding to the intersection as the maximum frequency shift; wherein, the slope of the straight line is correlated to the noise level; and the noise is estimated by a method similar to the percentage method, i.e., averaging the spectral lines adjacent to a cut-off frequency.
In an article “Comparison of the performance of three maximum Doppler frequency estimators coupled with different spectral estimation methods”, Ultrasound in Med. & Biol., Vol. 20, No. 7, pp 629-638, (1994), K. Marasek et al made an improvement on the mixed method, and proposed a geometric method comprising calculating the spectrum integral curve of each of the spectral line first, then analyzing each of the integral spectral lines. This method differs from the mixed method in the method for determining the maximum frequency shift point, in which a straight line is designed, the distances from each of the points on the integral spectral line to the straight line is calculated, and the frequency corresponding to the point with the minimum distance is regarded as the maximum frequency shift.
In an article “The performance of three maximum frequency envelope detection algorithms for Doppler signals”, J. vasc. Invest, 1:126-134 (1995), R. Moraes et al made an improvement on the geometric method, in which a straight line passing through the origin is designed, and the frequency corresponding to the point with maximum vertical distance between the integral curve and the straight line is regarded as the maximum frequency shift. In order to prevent inestimable errors from occurring when the signal is weak, and experimental threshold is additionally used in this method. If the signals are weaker than the threshold, the detection of maximum frequency shift is not performed, and the maximum frequency shift is directly set to the wall filtering cut-off frequency.
In the technical solution disclosed in U.S. Pat. No. 5,287,753, Routh et al introduced an adaptive threshold envelope detection method. The basic idea thereof comprises determining the maximum frequency shift by comparing the Doppler power spectral intensity with a set threshold. In this method, the threshold can be adjusted adaptively based on the SNR of each of the cardiac cycles. This method differs from the above mentioned methods in that the method assumes the noise level and mean SNR to be relatively stable in a cardiac cycle, therefore, the mean noise level and mean SNR are determined with respect to a cardiac cycle. However, the thresholds of the previous methods acutely vary with each of the spectral lines, and thus estimated threshold levels will also vary acutely with each of the spectral lines.
The main disadvantage of the above-mentioned techniques is that: the mixed method is based on the hypothesis that the maximum frequency shift in the frequency spectrum diagram shall be less than the cut-off frequency. The mean estimated noise in the boundary portion of the spectrum diagram is significative only under this hypothesis; thus this method will be invalidated when the cut-off frequency is less than the maximum frequency shift. In addition, since the result of maximum frequency shift estimation is greatly affected by the actual estimated noise and the integral curve, when the noise threshold is set to a small value, the maximum frequency shift estimation will be relatively small, and when the noise threshold is set to a big value, the maximum frequency shift estimation will be relatively big; thus this method is greatly affected by the actual SNR and accurate estimation of noise.
In the geometric method, since the maximum frequency shift estimation is a biased estimation, the effectiveness of this method will be lowered when the cut-off frequency is less than the maximum frequency shift. This method does not estimate the noise level directly, but it is severely affected by the integral curve; when the SNR is relatively low, the integral curve fails to guarantee a relatively rapid variation when the frequency is low, and a slow variation when the frequency is high, the error of the maximum frequency estimated then will be relatively high.
In the improved geometric method, since the maximum frequency estimation is a biased estimation that is always less than the actual value, similar to the geometric method, the estimation error will be greater when the maximum frequency shift approaches the cut-off frequency. Although, in this method, a threshold is added to lower misjudgement when the signal is weak, but the selection of the threshold depends on experiments; the differential of selections is relatively great under different conditions of SNR. Similar to the geometric method, when the signal energy is relatively high and SNR is relatively low, the integral curve fails to guarantee that it varies relatively rapid when the frequency is low, and relatively slow when the frequency is high, the error of the maximum frequency estimated then will be relatively high.
In the adaptive threshold method, since the noise threshold is obtained by averaging the SNRs, and the SNRs are relatively severely affected by the signal intensity, the maximum frequency shift estimation will be affected by an incorrect estimation of the noise level. Further more, since the threshold is in relation to the result of the previous cardiac cycle envelope detection, when the Doppler spectrum diagram is relatively stable, the adjustment of the adaptive threshold can not guarantee that the threshold will converge to a stable value. A simple threshold determination method is used in the determination of the envelope, the robustness of the envelope detection will also be impaired by the relatively great undulation of the power spectrum calculated by FFT.